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Abstract. We obtained exact heat capacities of the quantum compass model on the square 
L X L clusters with L = 2, 3, 4, 5 using Kernel Polynomial Method and compare them with heat 
capacity of a large compass ladder. Intersite correlations found in the ground state for these 
systems demonstrate that the quantum compass model differs from its classical version. 



Quantum compass model originating from the orbital (Kugel-Khomskii) superexchange in 
transition metal oxides has been recently studied both using analytical [T] and numerical [21 [3] 
approach. In spite of remarkable progress in numerical methods for two-dimensional (2D) Ising- 
like models [HE], exact solutions are necessary to investigate the structure of excited states. The 
one— dimensional (ID) compass model is exactly solvable by mapping to quantum Ising model 
[H] and exhibits interesting properties while approaching to the quantum critical point at zero 
temperature [7j. In addition, its ladder version, first considered by Dougot et al. [IJ is solvable 
in a similar way and its partition function [H] can be obtained exactly in case of a large (but 
finite) system. There is no exact solution for the 2D compass model but the latest Monte Carlo 
data [3J prove that the model exhibits a phase transition at finite temperature both in quantum 
and classical version with symmetry breaking between x and z part of the Hamiltonian. This 
paper suggests a scenario for a phase transition with increasing cluster size by the behavior of 
the specific heat obtained via Kernel Polynomial Method (KPM) 

The Hamiltonian of the quantum compass model on a square L x L lattice is given by 

L 

n{a) = (-l)^J E {(1 - «)<-^^?+i,«; + a<^<»+i} , (1) 

i,w=l 

where o"^'^ are the z and x Pauli matrices for site where i {w) is a vertical (horizontal) 

index, and we implement periodic boundary conditions. The sign factor (—1)^ is introduced 
to provide comparable ground state properties for odd and even systems. Parameter q E [0, 1] 
makes this model interpolate between the situation when we have L independent Ising chains 
interacting with x and z components of spins for a = and a = 1 respectively. The case we 
have already discussed is the compass model on a ladder [8]; this can be included into present 
discussion by restricting the range of the summation over i to the value 2 in the Hamiltonian. 

The latter case is especially interesting, as we can obtain exact ground state for any size L of 
the system. The quantities referring to the full spectrum, like the density of states, heat capacity 
and thermodynamic correlation functions can be determined for L sufficiently big (like L = 52) 



to be representative for the thermodynamic hmit. Solution for a ladder system is based on the 
construction of invariant subspaces which are related to the symmetry operators erf ^crfu,- It 
brings us to a purely ID Hamiltonian describing an Ising chain in transverse field but depending 
on in which subspace we are — some of the crfaf_^i bonds are missing. An analogous construction 
is possible for a square lattice but in this case simplifications are rather modest; we cannot find 
an exact solution anyway. This is not the case for a finite system - exact diagonalization (ED) 
methods can be applied and symmetries can reduce the Hilbert space considerably. 

Ground state energies and energy gap of the model given by ([T]) has been already calculated 
for different values of a and for L G [2, 5] using ED and for higher L using Green's function 
Monte Carlo method [2]). Our approach will be based on KPM [U] which will let us calculate 
the densities of states and the partition functions for square lattices of the sizes up to L = 5. 
We start by applying Lanczos algorithm to determine spectrum width which is needed for KPM 
calculations. The resulting few lowest energies that we get from the Lanczos recursion can be 
compared with the density of states to check if the KPM results are correct. One should be 
aware that the spectra of odd systems are qualitatively different from those of even ones. For 
the even systems operator S defined as = nf«)=i ^{^ ~ {~^y~^^}'^i anticommutes with the 
Hamiltonian ([T|). This means that for every eigenvector \v) satisfying Tl(a)\v) = E{a)\v) we 
have another eigenvector \w) = S\v) that satisfies TC{a)\'w) = —E{a)\w) . This proves that even 
values of L spectrum of TC{a) is symmetric around zero but for odd L's this does not hold; S 
no longer anticommutes with the Hamiltonian. To obtain symmetric spectrum we would have 
to impose open boundary conditions. 

We would like to highlight that KPM calculation for 5x5 lattice (2^5 - dimensional Hilbert 
space) would be impossible without using the symmetry operators. These operators are usually 
called Pi and Qw [1\ and are defined for i, tt; = 1, 2, • • • , L as 

L L 

P^=l[ Q«, = n (2) 

w=l i=l 

One can easily check that although both operators commute with the Hamiltonian, Pi and Quu 
anticommute. Thus, we cannot find a common eigenbasis for P's, Q's and the Hamiltonian; we 
should find a different set of operators. A good choice is to take all Pi with i = 1,2, . . . , L and 
Rw = QwQw+i with w = 1,2, . . . , L — 1. This gives us (L — 1) x L commuting symmetries. Let's 
denote their eigenvalues as pi and taking pseudospin values ±1. Each choice ol pi,p2, ■ ■ ■ ,Pl 
and ri,r2, • • • ,rL^i defines an invariant subspace in which the Hamiltonian can be written in 
terms of (L — 1) x (L — 1) new spin operators and (L — 1) x L parameters {pi} and {r^}. This 
statement can be proved by giving the explicit form of the spin transformation. 

The main benefit for us is that after the transformation the Hamiltonian of 5 x 5 compass 
model (a = 1/2) turns into 2^ spin models, each one on 4 x 4 lattice. In fact, the number of 
different models is much lower than 2^; most of resulting Hamiltonians differ only by a similarity 
transformation. We can check this using Lanczos algorithm; if the two lowest energy levels from 
two subspaces are the same then it is reasonable to assume that these spectra are identical and 
the Hamiltonians are the same. Finally, we find out that only 10 out of 512 Hamiltonians are 
different; their two lowest energies and degeneracies ars given in table [H In fact, these energies 
are known with much higher precision (lO"*^) than that given in the table 1, and we also get quite 
good estimation for the highest energies. This gives us a starting point for KPM calculations. 

Kernel Polynomial Method is based on the expansion into the series of Chebyshev polynomials 
[9]. Chebyshev polynomial of the n-th degree is defined as Tn{x) = cos[n arccos x] where 
X G [—1, 1] and n is integer. Further on we are going to calculate r„ of the Hamiltonian so 
first we need to renormalize it so that its spectrum fits the interval [—1,1]. This can be done 
easily if we know the width of the spectrum. Our aim is to calculate the renormalized density of 




Figure 1. Heat capacities as functions of temperature for a = 2 the compass clusters Lx L 
of sizes L = 2,3 (panel (a), dashed-dotted and long-dashed lines) and L = 4, 5 (panel (b), dashed 
and solid lines) and for the compass ladder of 104 spins (black dotted line). 



states p{E) given by p{E) = (l/D) J2n=o ^(E — En), where the sum is over eigenstates of 7i{a) 
and D is the dimension of the Hilbert space. The moments pn of the expansion of p{E) in basis 
of Chebyshev polynomials can be expressed by 

Pn = Tn{E)p{E)dE = i T,{Tn{rL)] . (3) 

Trace can be efficiently estimated using stochastic approximation: 

Tr{r4W)}«if:(r|r„(W)|r), (4) 

r=l 

where \r) (r = 1,2, ... ,R) are randomly picked complex vectors with components Xr,k (^ = 
1,2,..., D) satisfying {xr,k) = 0' {Xr,kXr',i) = 0, {Xr,kXr\i) = K,r'h,i (average is taken over 
the probability distribution). This approximation converges very rapidly to the true value 
of the trace, especially for large D. Action of the Tn{7i) operator on a vector |r) can be 
determined recursively using the following relation between Chebyshev polynomials: Tn(7i)\r) = 
{2'HTn-ii'H) - Tn-2{'H)}\r) . We can also use the relation 2Tm{x)Tn{x) = Tni+n{x) + Tm-n{x) 
to get moments p2n from the polynomials of the degree n. Finally the required function, 

P{E) « _ I qqPo + 2 ^ gnPnTn{E) I , (5) 



Table 1. Ground state energy £^0 and first excited state energy Ei and their degeneracies d for 
10 nonequivalent subspaces of the 5x5 compass model (1) at a = 1/2. 
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Figure 2. Correlations {(Tf^af_^_i ^) in 
the ground state calculated for sizes L = 
2, 3, 4, 5 (dashed-dotted, long-dashed, dashed 
and solid line) and for infinite ladder (black 
dotted line) as functions of a. Step 
function shows a discontinuous transition for 
a classical compass model at a = 1/2. 



can be reconstructed from the known moments, where gn coefficients come from the integral 
kernel we use for better convergence. Here we use Jackson kernel. Choosing the arguments 
of p{E) as being equal to = cos[(2A; — l)7r/2A^'] (k = 1,2, .. . ,N') we can change the last 
formula into a cosine Fourier series and use Fast Fourier Transform algorithms to obtain rapidly 
p{Ek). This point is crucial when and N' are large, which is the case here; our choice will be 
A'^ = 20000 and A^' = 2A'^. In such a way we can get the density of states for 4x4 and 5x5 
systems. In the latter case we obtain 10 energy spectra for 10 nonequivalent subspaces — these 
can be summed up with proper degeneracy factors (see table [1]) to get the final density of states 
p{E) and next the partition function via rescaling and numerical integration. 

The heat capacity Cy for the compass L x L clusters behaves differently from that for 
a compass ladder, see figure [H The main difference is vanishing of the low-T peak when the 
system's size increases. This correspond to vanishing of the low-energy modes which is consistent 
with presence of the ordered phase for finite T in the thermodynamic limit. In contrast, heat 
capacity of the ladder indicates robust low-energy excitations and dense excitation spectrum at 
higher energies causing a broad peak in Cy. 

In figure [2] we compare nearest neighbor correlations as functions of a for finite clusters, 
infinite compass ladder and classical compass model on a square lattice. Curves for finite clusters 
converge to certain final functions which is something intermediate between classical case and 
quantum ladder. This result shows that even in large L limit the 2D compass model preserves 
quantum correction to a classical behavior even though it chooses ordering in one direction [3]. 

In summary, we have shown that exact heat capacities of square compass Lx L clusters could 
be obtained by implementing the symmetries up to L = 5. The behavior of the low-T peak 
in the heat capacity indicates that the gap in the spectrum decreases with increasing L. This 
agrees with the numerical results obtained before by numerical approach [2j. 
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